Nonzero spontaneous electric polarization in metals: novel predictive methods and applications

Ferroelectricity in metals has advanced since the initial discovery of nonmagnetic ferroelectric-like metal LiOsO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}3, anchored in the Anderson and Blount prediction. However, evaluating the spontaneous electric polarization (SEP) of this metal has been hindered by experimental and theoretical obstacles. The experimental challenge arises from difficulties in switching polarization using an external electric field, while the theoretical limitation lies in existing methods applicable only to nonmetals. Zabalo and Stengel (Phys Rev Lett 126:127601, 2021, 10.1103/PhysRevLett.126.127601) addressed the experimental obstacle by proposing flexoelectricity as an alternative for practical polarization switching in LiOsO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}3, which requires a critical bending radius similar to BaTiO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}3. In this study, we focus on resolving the theoretical obstacle by modifying the Berry phase and Wannier functions approaches within density functional theory plus modern theory of polarization. By employing these modifications, we calculate the SEP of LiOsO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}3, comparable to the polarization of BaTiO\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}3. We validate our predictions using various ways. This study confirms the coexistence of ferroelectricity and metallicity in this new class of ferroelectric-like metals. Moreover, by addressing the theoretical limitation and providing new insights into polarization properties, our study complements the experimental flexoelectricity proposal and opens avenues for further exploration and manipulation of polarization characteristics. The developed approaches, incorporating modified Berry phase and Wannier function techniques, offer promising opportunities for studying and designing novel materials, including bio- and nano-ferroelectric-like metals. This study contributes to the advancement of ferroelectricity in metals and provides a foundation for future research in this exciting field.


Band classification and introduction to polarization methods for LiOsO 3
Both modified polarization methods, mBp and mWf, begin with band classification as their initial step.We will first cover this common step before individually discussing the subsequent steps for each method.To initiate our methodology, we will classify the bands of LiOsO 3 as calculated using the PBE-GGA functional for its rhombohedral structure, illustrated in Fig. 1.

Classification of bands
Classification of the valence and conduction bands constitutes the zeroth step of our modified methods of polarization.
In the zeroth step, we classify the bands into three classes, as shown in Fig. 1.These classes are labeled by I, I * and II in Fig. 1.These labels are dual-purpose and, in addition to their roles in classifying the bands, they can be also used to distinguish three different energy intervals indicated in Fig. 1.The class I includes valence bands which are crossing the Fermi level and conduction bands.The highest energy limit of class I * is the Fermi energy, while its lowest energy limit possesses in common with that of class I. Therefore, the class I * is a subclass of class I.Only some of the valence bands in classes I and I * are full while their other remaining valence bands are partially filled and hence they can only partially contribute into the electronic part of the Berry phase, as expressed in Eq. ( 12) of the SMs.The conduction bands in class I are empty and hence their contributions into the electronic part of the EP are zero.The class II includes valence bands which are all fully occupied.The valence bands included in class II are well separated by an energy interval of ≈ 0.4 eV from the bands included in class I, see Fig. 1.
After the classification, we introduce our polarization methodologies: the mBp and mWf methods.We begin with the six-step mBp polarization method.For a more detailed view, refer to Sect. 3 of the SMs.

Mean-field-like mBp method of polarization
To overcome the obstructs, as discussed in Sect. 2 of the SMs, and calculate the SEP of a FE-L metal, we modify the conventional Berry phase method, beginning after the zeroth step.The critical points of the proposed methods is concisely provided in the next six steps.For a comprehensive discussion, refer to Sect.3.1 of the SMs.
In the first step, we obtain the ionic part of the Berry phase, ϕ ( ) ion,µ , for the structure along polarization direction µ by applying the standard Berry phase approach over the density functional theory (DFT) results.
Here, , which can be 0 or 1, denotes the initial non-polar CS R 3 c structure and the final polar NCS R3c phase, which has a lower symmetry than R 3 c; i.e., ( = 0) ≡ R 3c and ( = 1) ≡ R3c.
In the second step, using the standard Berry phase approach, we obtain the electronic part of the Berry phase in the µ direction for all bands of class II at every perpendicular wave vector k ⊥ .This is illustrated in Fig. 2 for the structure , represented as ϕ ( ),(II) el,µ (k ⊥ ) .We then adjust the electronic Berry phase of class II to be independent of k ⊥ .To achieve this, we draw an analogy to the standard Berry phase method.As depicted in Fig. 2, we take an average either over the perpendicular area A ⊥ or over the discrete points N k ⊥ of the 2D k ⊥ -point samples.This allows us to determine ϕel, µ ( ),(II) as follows: Figure 1.(a) Band structure of the polar NCS R3c phase of the LiOsO 3 compound constructed by the rhombohedral (trigonal) crystal system calculated using the PBE-GGA 69 .(b) The corresponding first Brillouin zone including the symmetrical points and selected paths, as plotted by XCrySDen package 70 .
In the third step, we determine the electronic part of the Berry phase for all bands of class I, represented as ϕ ( ),(I) el,µ (k ⊥ ) .We make a temporary assumption that all bands of this class are fully occupied.This assumption will be modified in the fourth step.At this stage, we delay the averaging process for ϕ ( ),(I) el,µ (k ⊥ ) until the end of the fourth step.This is when the modification occurs.Hence, for now, it remains dependent on k ⊥ .
In the fourth step, we adjust ϕ ( ),(I) el,µ (k ⊥ ) by taking into account the correct occupation numbers for the class I bands, denoted as n ( ) n (k ⊥ ) , as follows: where φ stands for the modified ϕ , and M I is the number of bands of class I, as well as N ( ) k BZ is the total number of k-points generated in the full mesh of the first Brillouin zone for structure , see Fig. 2. Now, it is time to make φ ( ),(I * ) el,µ (k ⊥ ) independent of k ⊥ using the same averaging procedure outlined in Eq. (2) to derive φ ( ),(I * ) el,µ .In the fifth step, we combine the electronic Berry phases calculated in steps 2 and 4 to obtain the modified total electronic Berry phase for structure in direction µ .This is given by: φ ( ) el,µ = ϕ ( ),(II) el,µ + φ ( ),(I * ) el,µ .To find the total Berry phase for structure in direction µ , denoted as φ ( ) µ , we add the electronic Berry phase φ ( ) el,µ to the ionic Berry phase for structure in direction µ that was calculated and stored in step 1, represented as ϕ ( ) ion,µ .That is, φ ( ) µ = 2φ ion,µ , where the factor of 2 accounts for the spin degeneracy in non-spin-polarized systems.For spin-polarized systems, the expression is modified as: φ ( ) µ = ϕ , where (↑) and (↓) denote spins up and down, respectively.
In the sixth step, we substitute ϕ ( ) µ into the following equation: (1) ϕ ( ),(II) el,µ el,µ (k ⊥ ), (2) φ ( ),(I * ) el,µ where e is the electron charge and � ( ) is the unit cell volume of the structure .In Eq. (3), R ( ) µ is the length of the lattice vector in real space for the structure , viz.R ( ) = 3 µ=1 R ( ) µ ê( ) µ , where R ( ) µ ê( ) µ (R ( ) µ ) is the primitive vector (lattice constant) of structure along µ in the direction of the unit vector ê( ) µ .Then, after multiplying both sides of Eq. ( 3) by êµ and taking the summation over µ on both sides of the resultant equation, we obtain the polarization vector for structure as: The procedure discussed above, from step 1 to this stage of step 6, is performed for both structures: one for structure = 0 and the other for structure = 1 .In this manner, we obtain the electric polarization vectors P ( =0) for the structure = 0 and P ( =1) for the structure = 1.

mWf method of polarization
The conventional Wannier functions scheme is not apt for FE-LMs as it presumes all bands to be fully filled, a condition not met by FE-LMs.Consequently, we employ the partly occupied maximally localized Wannier functions methodology introduced by Thygesen et al. 71,72 , further adapting the occupation numbers beginning after the zeroth step.
This method has shown that the degree of localization of Wannier functions can be optimized with a specific number of unoccupied orbitals 71,73 .Andrinopoulos et al. 74 enhanced DFT's van der Waals energy contributions using partly occupied MLWFs, incorporating anti-bonding states.For metals, considering only occupied states can result in poorly localized Wannier functions.However, adding unoccupied conduction states can significantly improve localization 65,71,73,74 .Although the maximally localized Wannier functions have been extended for partly occupied states 71,72,74 , their application for predicting electric polarization in FE-LMs is unexplored.Occupancy in polarization formulas is often assumed to be 2 for each Wannier center, but metals with partial occupation can have values less than 2. We find it pertinent to provide a succinct overview of this method for SEP calculations in FE-LMs, incorporating modifications to tackle these challenges.For an exhaustive discussion, readers are directed to Sect.3.2 of the SMs.
In our research, we used the approach delineated by Thygesen et al. 71 from 2005.This methodology is wave vector-independent, paving the way for its straightforward implementation in both isolated and periodic systems.On the other hand, an earlier method presented by Souza et al. 61 in 2001 uses a 'disentangling procedure' .This technique, while older, zeroes in on specific bands, aiming to reduce variations in the character of Bloch states across the Brillouin zone, with an emphasis on revealing pertinent unoccupied states.Although these methods vary, their shared objective revolves around the creation of more localized Wannier functions by examining the conduction and valence bands proximate to the Fermi level.In contexts where band groupings are not evident or when specific computational attributes of bands near the Fermi level are required (as in the calculation of the anomalous Hall conductivity via Wannier interpolation 75 ), the method by Souza et al. 61 might be the go-to choice.A testament to its utility is the work of Wang et al. 75 , where they applied a post-processing step to Bloch states near the Fermi level.Leveraging the Souza et al. technique 61 , they mapped these states onto localized Wannier functions, enabling the computation of the anomalous Hall conductivity.This approach facilitated precise Berry curvature interpolation for any selected k-point, striking a balance between computational efficacy and precision.However, in scenarios echoing our situation, where distinct band groups are discernible (as depicted in Fig. 1), the partly occupied approach by Thygesen et al. 71 is potentially more beneficial.Hence, for the purposes of our study, we deem the newer partly occupied method 71 to be aptly suited.
The total polarization vector P ( ) for structure can be expressed in terms of its electronic P ( ) el and ionic P ( ) ion contributions within the Wannier functions framework as: where the first (second) term can be interpreted as ionic (electronic) polarization per unit volume � ( ) of struc- ture originated from N ions ( J Wannier centers) each with positive (negative) charges of +eZ is the occupancy of Wannier center n and J is the number of Wannier centers in structure .
In the ionic polarization, r ( ) s represents the classical position of ion s in structure .On the other hand, in the electronic polarization, r W ( ) n,R signifies the expectation value of the position of the Wannier center n in (3) (5) , structure .It's essential to understand that r Wn, R ( ) does not correspond to the position of a classical particle and should be interpreted within a quantum mechanical framework.
In order to calculate P ( ) , we first decompose the electronic polarization based on the number of composite bands.For LiOsO 3 , there are two distinct composite bands: • The isolated class of bands II, which includes only deep-lying, fully occupied valence bands.
• The isolated class of bands I, which comprises shallow-lying, fully and partially occupied valence bands, as well as low-lying empty conduction bands.
Thus, we decompose P el at lattice vector R in structure .The factor 2 inside the first summation accounts for fully occupied Wannier centers, whereas n ( ) W n,R inside the second summation accounts for both fully and partially occupied Wannier centers.We have J = J (II) + J (I) .The term J (I) , as the upper limit of the second sum in Eq. ( 6), refers to the Wannier center of class I.However, the occupation numbers n W ( ) n,R used in the second term of Eq. ( 6) are so determined subsequently in step three that the second term itself refers to the polarization of class I * and yields P ( ),(I * ) el .To complete the calculation of P ( ) from Eqs. ( 5) and ( 6), we perform the following 5 steps.For a more detailed discussion on each step, readers are referred to Sect.3.2 of the SMs.
In the first step, we perform a regular self-consistent DFT calculation for the structure .Then, we restrict the energies to the energy interval II, as shown in Fig. 1.Now, we apply self-consistently the standard maximally localized Wannier functions procedure on the fully occupied valence composite bands of class II for structure , see Fig. 1.This procedure is performed over the Bloch states calculated by WIEN2k package 76,77 to obtain maximally localized Wannier functions and their centers of charges using Wannier90 code [60][61][62][63]65 and WIEN-2WANNIER interface 64 . At his stage, we use the first term of Eq. ( 5) to calculate P ( ) ion entirely while we use the first term of Eq. ( 6) to calculate partial electronic polarization P ( ),(II) el . Th latter electronic polarization is partial because it needs to be completed by including contributions of the Wannier centers of class I, i.e.P ( ),(I) el as the second term of Eq. ( 6).
In the second step, we restrict the energies to the energy interval I, see Fig. 1.Then, we apply self-consistently the generalized maximally localized Wannier functions procedure constructing partly occupied Wannier functions on the composite bands of class I including valence and conduction bands for structure , see Fig. 1.By this way, we calculate the positions of the Wannier centers of the structure as for Wannier centers n = 1 to J (I) .At this stage, we cannot calculate the remaining electronic polarization by the following conventional equation: because in this equation the occupation numbers of all the Wannier centers are assumed to be 2, while the class I contains partially occupied Wannier centers.We cannot also use our generalized formula expressed as the second term of Eq. ( 6), because n W ( ) n,R still are unknown.Therefore, the main task of the next step is to determine the unknown occupation numbers n W ( ) for the Wannier centers of class I.
In the third step, we determine n W ( ) n,R so that the polarization calculated in the following fourth step gives the polarization P ( ),(I * ) el which is related to the desired class I * .To this end, we first individually project the density of states (DOS) on each of the maximally localized Wannier centers for n = 1 to J (I) , see Fig. 3a.For instance, the projected DOS, as shown in Fig. 3a, is obtained by projecting the calculated DOS on one of the maximally localized Wannier centers which is shown in Fig. 3b.Then, we integrate each of the projected DOSs up to the Fermi level.By this, we obtain individually the areas under each of the DOSs projected on the maximally localized Wannier centers up to the Fermi level, e.g.see the filled area under the DOS shown in Fig. 3a n,R into the second term of Eq. ( 6), P el for structure using the generalized Eq. ( 6).Let us close this step by indicating a practical note.To this end, let us consider a system whose its polarization direction is pointed along only a single direction.Such a system resemblances the case under study whose polarization direction is oriented along the c axis of the hexagonal supercell.In this case, it is enough to consider only n W ( ) In the fifth step, based on Eq. ( 5), we add the ionic part of polarization P ( ) ion for structure , as obtained in the first step, to the electronic part of polarization P ( ) el for structure , as obtained in the fourth step.By this, we obtain the total electric polarization P ( ) for structure .

Merits, limits, and management of mBp and mWf methods of polarization
Here, let us assess the strengths and limitations of our mBp and mWf methods of polarization.While most of computational methods inherently have distinct advantages due to their foundational principles and algorithms, they can also face certain challenges or constraints.Rooted in specific theoretical frameworks, our mBp and mWf methods have been developed to offer particular strengths designed for certain applications.Nevertheless, they are not without limitations.In the subsequent sections, we detail and address these constraints.Our goal is to provide a thorough understanding of these methods' scope and to highlight situations where they demonstrate optimal effectiveness.mBp method Our mBp methodology represents a detailed evolution in the filed of polarization calculations.At its core, the mBp approach utilizes mean field-like calculations, with a focused attention on determining the occupation weight, predominantly for k ⊥ .The merit of this approach is most evident when applied to materials like LiOsO 3 , which exhibit distinct isolated band structures.
One of the notable features of our method is its adaptability to the material's intrinsic electronic configuration, making it versatile for various materials, including metals and insulators.We observe that the SEPs calculated by mBp (PBE-GGA) and Bp (PBE-GGA+U) for LiOsO 3 are close to each other for the material in question, showing the limited impact of bands crossing the Fermi level on the resultant polarization in this specific instance.This observation validates the applicability and reliability of our methodology for systems that exhibit the electronic behavior similar to LiOsO 3 .www.nature.com/scientificreports/However, it is essential to understand the boundaries of any method's universality.In scenarios devoid of such isolated band structures, our mBp might need further refinements.One promising direction involves integrating the Wannier bases approach to compute the integral tied to the Berry phase, drawing inspiration from Wang et al. 75 .For cases influenced by d-orbitals around the Fermi level, the Hubbard model, as deduced from our prior research 78 and validated in this study, emerges as a powerful tool, enabling us to simulate a gap and further analyze the SEP using the conventional mB without modifications.
To address the inherent challenges associated with the mBp method, particularly regarding the use of mean field-like calculations for occupation weight focused on k ⊥ , we incorporated a dense k-mesh in our post-pro- cessing calculations.Additionally, by comparing the occupation numbers for class I * bands-obtained by aggre- gating weights of k ⊥ -with those from self-consistent DFT calculations, we found noteworthy alignment.This compatibility reaffirms the reliability and soundness of our methodology.
Moreover, by emphasizing the change in polarization between CS and NCS phases, and consistent application of mean field-like calculations to both, we benefit from an inherent error compensation mechanism.
It is worth noting that in the conventional Bp method, the adiabatic condition is crucial for deriving Eq. ( 10) of the SMs.As our modification builds upon this equation, clarifying the physical rationale behind our adjustments becomes paramount.
Our mBp computational approach is designed to emulate conditions typical of insulating systems.As detailed in step 3 of "Mean-field-like mBp method of polarization" and further expanded in Sect.3.1 of the SMs, our calculations, particularly when using Eq.(10) (SM), align with the behavior of insulators.This choice ensures adherence to the adiabatic condition inherent in insulating systems.
To enhance the robustness of our method, we have incorporated refinements, such as considering weights at each k-point to adjust the previously calculated phase.These methodological tweaks aim to preserve an adiabatic-like behavior in our calculations, even if the system does not strictly abide by the adiabatic condition.
The underlying physical foundation of these adjustments can be understood as follows: in the adiabatic framework, external parameters like applied electric fields or strain, as the latter used in our case, vary slowly.This variation induces a change in the electronic polarization within a crystal due to adjustments in the selfconsistent Kohn-Sham potential.To capture this dynamic, we introduce a parameterization for the potential, which spans from 0 (initial potential) to 1 (final potential), covering 9 intermediate potential stages, as illustrated in Fig. SM6 of the SMs.
In summation, our mBp methodology, underpinned by rigorous scientific principles and augmented by integrated techniques, is a robust tool in polarization calculations.Its application across a spectrum of materials, judiciously taking into account its foundational strengths and suitable augmentations, promises reliable outcomes.

mWf method
Our mWf methodology, designed to address electronic structures like the one observed in our primary case, optimizes the post-processing approach, especially when there are distinct, isolated groups of bands, as visualized in Fig. 1.For such systems, our adaptation of the partially occupied approach by Thygesen et al. 71 is highly effective, achieving a consistent, reliable representation of electronic behaviors.
However, we acknowledge that in more intricate electronic landscapes, where there is no conspicuous energy window around the Fermi level, challenges can arise.In these instances, our method's intrinsic flexibility allows the incorporation of refined techniques, like the ' disentangling procedure' advocated by Souza et al.This procedure ensures that the smallest spread for subsequent Wannier functions is secured, bolstering the accuracy of our calculations.
Moreover, the electronic characteristics of our study material, prominently featuring d-orbitals around the Fermi level (elaborated in Figure SM2 and Sect.3.2 of the SMs), allow our mWf method to achieve optimal results with minimal iterations.It is paramount to note that while this property aids our specific case, there could be materials where orbital hybridizations are more pronounced.For such cases, our approach remains versatile: by carefully selecting various hybrid orbitals as initial estimates, we ensure the generation of highly localized Wannier functions, even in the face of significant hybridization challenges.
In summary, we have designed our mWf method to be both adaptive and resilient, able to cater to a diverse range of electronic structures while maintaining a high degree of precision.By recognizing potential challenges and proactively integrating solutions into our methodology, we remain confident in the method's applicability and accuracy across varied electronic landscapes.

SEP direction in LiOsO 3 FE-LM
First, it would be sensible to differentiate the non-polar phase of LiOsO 3 FE-LM from its polar phase and subsequently determine qualitatively the direction of the spontaneous electric polarization of the system under consideration.For the sake of conciseness, discussions on the CS R 3 c and NCS R3c crystal structures of LiOsO 3 are detailed in Sect. 5 of the SMs, as illustrated in Fig. SM5.Notably, the R3c structure can be characterized as the polar phase, denoted by " = 1 " representing the final structure of the ferroelectric-like phase transition, while the R 3 c structure serves as the non-polar phase, indicated by " = 0 ", representing the initial structure of the transition.Based on the displacement directions identified between the polar and non-polar phases, the spontaneous polarization vector is anticipated to align along the c-axis.An in-depth exploration of this topic is presented in Sect.3.3 of the SMs, as shown in Fig. SM4.The spontaneous polarization orientation deduced in this section is further corroborated quantitatively in "SEP of LiOsO 3 : mBp approach of electric polarization" and "SEP of LiOsO 3 : mWf approach of electric polarization", where both the non-polar CS R 3 c and polar CS R3c phases are assessed using two distinct methodologies.In the modern theory of polarization [51][52][53][54][55][56][57][58][59] , the spontaneous polarization P is defined as the integrated current flow along the distortion direction during an adiabatic transition from a non-polar CS ( " = 0 ") structure to a polar NCS ( " = 1 ") structure, or in other words, �P := P ( =1) − P ( =0) modulo eR/� 52,55-57 , see Sect. 1.1 of SMs.
For the system under consideration, " = 0 " refers to the R 3 c structure, while " = 1 " corresponds to the R3c.The spontaneous polarization P = P (R3c) − P (R 3c) , which represents the change in polarization during the phase transition from CS to NCS, is a more significant physical quantity than the absolute polarizations of the CS and NCS structures.Despite this, to demonstrate the consistency between the results calculated by our mBp and mWf approaches of polarization, we will present the individual partial polarizations for both phases in addition to P. " Following the first step of the mBp approach, as discussed in "Mean-field-like mBp method of polarization" and Sect.3.1 of SMs, we have calculated the ionic parts of the polarizations, ϕ ( ) ion,µ , for both the non-polar CS R 3 c and polar NCS R3c phases taking the hexagonal supercells containing 30 atoms into account.The results are presented in Table 1.Our numerical results, in agreement with the prediction discussed in "SEP direction in LiOsO 3 FE-LM", confirm that the directions of the polarizations are oriented along the c-axes of the hexagonal supercells.This means that ϕ ( ) ion,µ for µ = 1&2 are almost zeros and thereby negligible compared to ϕ ( ) ion,µ for µ = 3 .Hence, the µ index is known to be 3.Therefore, in Table 1, we omitted the known index µ , simplifying ϕ ( ) ion,µ to ϕ ( ) ion .The non-polar CS R 3 c (polar NCS R3c) is the initial (final) phase.The index for the initial (final) phase known to be 0 (1) refers to R 3 c (R3c), viz = 0 ≡ R 3c and = 1 ≡ R3c .Thus, for simplicity, the index is also omitted in Table 1, simplifying ϕ ( ) ion to ϕ ion .In the second step of the mBp approach, as discussed in "Mean-field-like mBp method of polarization" and Sect.3.1 of SMs, we have calculated the electronic part of the Berry phase for all the bands of class II, ϕ ( ),(II) el,µ , in both the non-polar CS R 3 c and polar NCS R3c phases.In analogy to the simplification made for the ionic part of Berry phase, ϕ ( ),(II) el,µ is also similarly simplified to ϕ (II) el , see Table 1.Applying the third and fourth steps with Eqs.(38) to (47), as discussed in "Mean-field-like mBp method of polarization" and Sect.3.1 of SMs, we have calculated the electronic part of the Berry phase for all bands of class I * , φ ( ),(I * ) el,µ , for both the non-polar CS R 3 c and polar NCS R3c phases.After omitting µ and indexes, as in the previous steps, the results are tabulated as φ (I * ) el in Table 1.Please, notice that φ differs from ϕ , see Eq. ( 38) and the notes after Eq. ( 39) of SMs, where φ is defined to be distinguished from ϕ.
Following the fifth step, we have first found the total electronic Berry phase, φ el , for both the phases indi- vidually by adding ϕ (II) el to φ el , see Table 1.We have then summed the total electronic Table 1.Calculated partial Berry phases and corresponding polarizations for the non-polar CS R 3 c and polar NCS R3c hexagonal supercells of LiOsO 3 , using the PBE-GGA with mBp and mWf schemes.The values for the electronic parts of the Berry phases ( φ el ) and polarizations ( P el ) were obtained by adding the contributions from class II ( ϕ (II) el and P (II) el ) and class I * ( φ el and P (I * ) el ).Similarly, total Berry phases ( φ ) and polarizations ( P ) were calculated by adding ionic and electronic contributions.Spontaneous polarizations ( P ) were derived by subtracting the polarizations of R3c and R 3 c structures as �P = P(R3c) − P(R 3c) .The results are given modulo eR/� and will be unwrapped in "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch" following the procedure proposed by Resta and Vanderbilt in Ref. 57 .Refer to Sect. 3 of the SMs for more detailed information.Berry phase φ el and the ionic Berry phase ϕ ion to find the total Berry phase φ = 2φ el + ϕ ion for both the phases individually, as presented in Table 1, where 2 shows the spin degeneracy.
In the sixth step, the total polarization P could be obtained by substituting the total Berry phase φ = 2φ el + ϕ ion into Eq. ( 2 of SMs.However, we here preferred to use the second way indicated in the sixth step to obtain not only the total polarization but also all the partial components of polarization.Therefore, following the second way of the sixth step, we have obtained the ionic polarization P ion by substitute the ionic Berry phase ϕ ion into Eq. (2)in Sect. 1.1 of SMs for both the phases individually, as tabulated in Table 1.Then, we have obtained the partial electronic polarizations P (II) el and P el , as well as the total electronic polarization P el by substitute the partial electronic Berry phases ϕ (II) el and φ el , as well as the total electronic Berry phase φ el into Eq.( 9) of SMs, receptively, for both the phases individually, as tabulated in Table 1.Consequently, we have obtained the total polarization as P = P ion + P el .
Eventually, using our mBp approach of polarization, we have calculated the spontaneous electric polarization �P(= P(R3c) − P(R 3c)) for the FE-LM LiOsO 3 , as reported in Table 1.
The results show that the ionic Berry phase for the CS (R 3 c) structure is 6.2832 rad , very close to an integer multiple of 2π rad (Table 1).Consequently, the ionic part negligibly contributes to the total polarization for the CS phase, as demonstrated by P ion = 0.0000 C/m 2 calculated for the R 3 c phase.
For the NCS phase, ϕ ion is 2.7396 rad , a value not equal to an integer multiple of 2π rad .Hence, according to Eq. ( 2) of SMs, this ionic Berry phase significantly affects the total polarization, leading to P ion = 0.3130 C/m 2 for the R3c phase.
According to Eqs. ( 2) and ( 9) of SMs, the ionic and electronic polarizations are obtained by multiplying ϕ ion 2π and φ el π by the quantum of polarization eR/� , respectively.ϕ (II) el is − 2.2402 rad for the polar R3c phase and 0.0000 rad for the non-polar R 3 c phase (Table 1).It indicates that the partial electronic Berry phase originating from the fully occupied deep-lying bands (class II) does not contribute to the polarization of the R 3 c phase, while it significantly contributes to the R3c phase.
Surprisingly, even in the non-polar CS phase, non-zero polarizations exist.For example, the φ el and φ el (= φ 1), resulting in In *"Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch", we examine the electronic Berry phases of LiNbO 3 , LiTaO 3 , BiFeO 3 , and LiOsO 3 and show that our calculated spontaneous polarizations are in agreement with existing experimental data and theoretical results.
However, despite the precise values obtained, we still have to account for the uncertainty rooted in the Berry phase theory of polarization, which defines polarization only modulo a quantum of polarization 57 .This suggests that polarization is a multivalued quantity.
In "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch", we address this uncertainty by calculating the polarization at several intermediate points along the transition path following the procedure by Resta and Vanderbilt 57 .This process allows us to select the best branch and to provide unwrapped results of the spontaneous polarization.

SEP of LiOsO 3 : mWf approach of electric polarization
As a part of the initial stage of the mWf approach, outlined in "mWf method of polarization" and Sect.3.2 of the SMs, we focus on the energy interval II.Within this context, we calculate the ionic component of polarization, represented as P ( ) ion , and the partial electronic polarization, denoted as P ( ),(II) el .We carry out these calculations with a focus on two distinct phases: the non-polar centrosymmetric (CS) R 3 c phase, and the polar non-cen- trosymmetric (NCS) R3c phase, taking into account the hexagonal supercells.The results are presented as P ion and P (II) el in Table 1, where the known indexes are removed for simplicity for both the CS R 3 c and NCS R3c phases.Similarly, the indexes µ have been eliminated because they consistently yield a value of 3 for both the CS R 3 c and NCS R3c phases.Our computational findings from the mWf method reveal that the x and y components of both P ( ) ion and P ( ),(II) el are essentially zero, meaning they are remarkably close to integer multiples of eR/� .Our results from the mWf method indicate that the x and y components of P ( ) ion and P ( ),(II) el are negligible compared to their respective z components in both the centrosymmetric (CS) R 3 c and non-centrosymmetric (NCS) R3c hexagonal supercells.This observation aligns with the numerical predictions of the mBp scheme (as discussed in "Mean-field-like mBp method of polarization" and Sect.3.1 of the SMs) and the theoretical predictions discussed in "SEP direction in LiOsO 3 FE-LM".This evidence substantiates that the polarization vectors P el for the CS R 3 c phase using the mWf method, align perfectly with their respective values calculated by the mBp scheme, as shown in Table 1.Moreover, the mWf-calculated value of the partial electronic polarization P ( ),(II) el for the CS R 3 c phase, measured at 0.2076 C/m 2 , is in close agreement with the mBp-calculated value of 0.2060 C/m 2 for the same phase (see Table 1).
Utilizing the second step of the mWf method, as detailed in "mWf method of polarization" and Sect.SM2 of SMs are calculated for the rhombohedral unit cells, the DOSs calculated for the hexagonal supercells, not presented here, show approximately similar behaviors.In the hexagonal supercell, there are six Os 5+ ions.Each Os 5+ ion has a nonmagnetic 5d 3 ground state, leading to almost 3 valence d-electrons per Os 5+ ion 17,20,26 .Our results, in agreement with Refs. 17,20,26, show that the metallic state of LiOsO 3 mainly originates from the d-orbital of the Os 5+ ions, see Fig. calculated below confirm that the x and y components of the Wannier centers corresponding to the region I do not contribute to the polarizations in both the CS R 3 c and NCS R3c hexagonal supercells.Thus, only the z components of the Wannier centers z W n,R for n = 1 to J (I) = 18 are tabulated here in Table 2 for both the CS R 3 c and NCS R3c phases.
Following the third step, we have determined n W ( ) for both the CS R 3 c and NCS R3c phases.To do this, we have first projected the total Wannier DOS on each of the 18 maximally localized Wannier centers for n = 1 to J (I) individually.Then, we have integrated each of the projected Wannier DOSs up to the Fermi level one by one.By integrating up to the Fermi level, we have changed the working class from the undesired I to the desired I * .The areas under the projected Wannier DOSs calculated up to the Fermi are tabulated in Table 2 as the occupation numbers n W n,R , after removing the known indexes, for n = 1 to J (I) = 18 .We have examined the correctness of the occupation numbers by summing on n W n,R over all the Wannier centers for both the CS R 3 c and NCS R3c phases individually.The examination, as also presented in Table 2, validates that  to the conduction states.The element A is repeated twice while the element B is repeated once in the AAB pattern.The amount of the occupation numbers can depend on the number of bands that crosses the Fermi level and the ratio of the number of conduction states added to the total number of valence and conduction states.
For the class I, there are 6 bands that cross the Fermi level, and 6 valence bands, as well as 6 conduction bands.The 12 Wannier centers related to the 6 valence bands and the 6 conduction bands are closer to the positions of the Os 5+ ions while the remaining 6 Wannier centers related to the 6 bands that cross the Fermi level are farther from the positions of the Os 5+ ions.The occupation number of the 12 Wannier centers which are closer to the ionic positions is larger than that of the remaining 6 Wannier centers which are farther from the ionic positions.The larger (smaller) occupation number is A = 1.09 ( B = 0.82 ).This results in the AAB pattern.Following the fourth step, we have determined the partial electronic polarizations P ( ),(I * ) el,µ for = 0 and 1 by substituting the multiplications of n W ( ) into the second term of Eq. ( 52) of SMs using the results r W ( ) tabulated in Table 2.The resultant partial electronic polarizations are given for both of the phases in Table 1 as el , where the known indexes and µ are removed.We have checked that the x and y components of this partial electronic polarization are very close to integer multiples of eR/� leading to vanished polarizations along x and y directions for both the phases individually.The partial polarizations P (I * ) el are calculated by the mWf method to be 0.0706 C/m 2 for the NCS R3c and − 0.3584 C/m 2 for the CS R 3 c which are close to the correspond- ing partial polarizations P (I * ) el calculated by the mBp method, i.e. 0.0746 C/m 2 for the NCS R3c and -0.3519 C/ m 2 for the CS R 3 c.This shows that both the mBp and mWf approaches yielding consistent results can be con- sidered as two different reliable methods to predict polarization corresponding to the entangled bands of class I * .Then, we have obtained the electronic polarizations P el for both of the phases using the mWf method expressed in the generalized Eq. ( 52) of SMs by adding P el , as tabulated in Table 1.The results show that the electronic polarizations P el calculated by the mWf method for both of the phases are in agreement with the corre- sponding polarizations calculated by mBp method, see Table 1.
Utilizing the fifth step of the mWf method, we have obtained the total electric polarizations P ( ) by substitut- ing P ( ) ion and P ( ) el , as tabulated in Table 1, into Eq.( 51) of SMs for = 0 and 1 .The total electric polarizations calculated by the mWf method are presented for both of the phases in Table 1 as P , where the known index has been removed.The results show that the total electric polarizations calculated by the mWf and mBf methods are consistent with each other, see Table 1.Eventually, we have obtained the spontaneous polarization P by �P = P ( =1) − P ( =0) according to the modern theory of polarization [51][52][53][54][55][56][57][58][59] .The spontaneous polarization P , as calculated by the mWf method, is 0.9496 C/m 2 which agrees with the value of 0.9455 C/m 2 calculated by the mBp method in "SEP of LiOsO 3 : mBp approach of electric polarization", see Table 1.This agreement authenticates that mWf and mBp are able to predict consistently spontaneous polarizations of FE-LMs.
Analogous to the spontaneous polarization of 0.9455 C/m 2 calculated by the mBp method in "SEP of LiOsO 3 : mBp approach of electric polarization", the value of 0.9496 C/m 2 obtained through the mWf method in this section is not considered the final result due to the quantum uncertainty problem.The phase freedom in the choice of the u nk , was shown to leave P el , invariant modulo eR/� 55 .The quantum uncertainty found in eR/� is reflected by the fact that the Wannier center position is defined only up to a lattice vector 79 .Therefore, the polarization can be considered as a multivalued quantity due to this uncertainty 79 .To overcome the quantum uncertainty problem of the mBp and mWf methods, the main task of the next section is devoted to counting the integer number of quanta involved in the polarizations calculated in "SEP of LiOsO 3 : mBp approach of electric polarization" and/or "SEP of LiOsO 3 : mWf approach of electric polarization".
In summary, the above discussion covers a multi-step computational method (the mWf approach) that deals with the calculation of ionic and partial electronic polarizations of the non-polar CS R 3 c and polar NCS R3c phases in certain hexagonal supercells.
In the first step, our calculations show that the x and y components of the polarizations are almost zero and therefore negligible in comparison to the z components.This implies that the polarization vectors are primarily aligned along the c-axis of the hexagonal supercells.Furthermore, these calculated values agree with prior calculations from the mBp scheme.
The second step involves calculating the positions of the Wannier centers, considering that there are 18 valence electrons predominantly originating from the Os atoms.It confirms that these 18 electrons are evenly distributed over the valence and low-lying conduction bands.Therefore, only the z components of the Wannier centers are considered significant and are tabulated.
The third step involves determining the occupation numbers of the Wannier centers by projecting the total Wannier DOS onto each center and then integrating up to the Fermi level.The occupation numbers are nearly equal to one, indicating that the 18 electrons are uniformly distributed over the 18 centers of the maximally localized Wannier functions.
The final step mentioned involves determining the partial electronic polarizations using the calculated occupation numbers and the positions of the Wannier centers from the previous steps.
Consequently, the mWf method accurately calculates the polarizations and verifies the orientation of these polarizations along the c-axis of the hexagonal supercells.It also calculates the positions and occupation numbers of Wannier centers.For more detailed information see Sect.In both the Wannier functions and Berry phase approaches of polarization, the spontaneous polarization P along an adiabatic path is a multivalued quantity that can be only well defined modulo a quantum of polarization eR/� 57 , where R is the lattice vector in the real space.In principle, there is such an uncertainty in polarization in both the Berry phase approach, as indicated in "Mean-field-like mBp method of polarization" and Sect.3.1 of SMs, and the Wannier approach, as indicated in "mWf method of polarization" and Sect.3.2 of SMs.In the Berry phase (Wannier) approach of polarization, a phase (Wannier center position) can be only well-defined modulo 2π ( R ).This implies that P can be defined uncertainly as P ( =1) − P ( =0) modulo eR/� 52,55-57 , which is a consequence of transnational symmetry 80 .The definition " �P := P ( =1) − P ( =0) (mod eR/�) " reads " P and P ( =1) − P ( =0) are congruent modulo eR/� ".This means that P and P ( =1) − P ( =0) can be different but equivalent in mod eR/� as they have the same remainder when divided by eR/� .In this definition, P is a factual quantity that can be observed and measured experimentally while P ( =1) − P ( =0) is a successor quantity proposed by the modern theory of polarization [51][52][53][54][55][56][57][58][59] that may not be necessarily equal to the factual quantity.In other words, computing P ( =1) − P ( =0) by the endpoints of the path only, may not always lead to the factual P .This is the case because there is no guarantee that the successor spontaneous polarization P ( =1) − P ( =0) is computed using the correct branch.If we only consider the endpoints of the path without verifying the branch's correctness, we might not obtain the accurate result 57 .Therefore, we have considered the uncertainty problem to uniquely obtain the spontaneous polarization of LiOsO 3 , as to be discussed subsequently.
Let us first more specifically clarify the problem.For the case under study, both of the polarizations P ( =0) and P ( =1) and consequently the spontaneous polarization P are oriented along the c axes of the hexagonal CS and NCS supercells, see "SEP direction in LiOsO 3 FE-LM", "SEP of LiOsO 3 : mBp approach of electric polarization", and "SEP of LiOsO 3 : mWf approach of electric polarization".Therefore, for this case, R employed in eR/� can be simplified as R = nc k so that |R| = R = nc , where n is an integer number and c ( k ) is the lattice constant (unit vector) along the Cartesian z axis.Hence, the above definition can be represented as  c,d).The auxiliary symbols ∧ and ∼ indicate that when wrapping and/or shifting are/is performed, if necessary, compared to the results presented in Table 1, see "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch" where the symbols are defined.All the Berry phases and as a result polarizations are calculated by the mBp scheme including non-spin-polarized PBE-GGA along the distortion path as functions of structure from " = 0 " to " = 1 " by step 0.1.The quantum of polarization and its number n are shown in (c,d).Our SEPs calculated by mWf, PBE-GGA+U with U = 0.2 and 2 eV are presented for comparison.The mBp, mWf, and empirical results presented in this figure are obtained for the metallic state of the NM LiOsO 3 .Our GGA+U results and the GGA+U result taken from Ref. 22 are calculated for the nonmetallic state of the G-AFM LiOsO 3 using the standard Berry phase method.The P Emp is extracted from Ref. 67 , generated using an empirical equation.The theoretical datum is taken from Ref. 22 .
�P := P ( =1) − P ( =0) + enc k/� or equivalently as �P k := P ( =1) k − P ( =0) k + enc k/� , where P = | P| , P ( =1) = |P ( =1) | , and P ( =0) = |P ( =0) | .By taking a dot product of the latter vector identity with the unit vector k , it can be simplified to its scalar form �P := P ( =1) − P ( =0) + enc/� .Therefore, the basic task to identify P uniquely is reduced to determine the integer number n for this case with polarization oriented along one- dimension only.We do it below by the procedure proposed in Ref. 57 .To this end, in addition to the starting structure = 0 " and end structure = 1 ", as the two endpoints of the adiabatic transition, we have constructed 9 intermediate structures = 0.1, 0.2, ..., 0.9 , as shown in Fig. C2 and discussed in details in Sect.5.4 of SMs.These intermediate structures are constructed using the freedoms of the structure = 1 ".The freedoms originate from the 5 internal parameters z1, and z2, as well as x3, y3, z3 existed in the potions of Li + , and Os + , as well as O 2− ions in the polar NCS R3c structure 81 , respectively, see Secs.5.1, 5.2, 5.3 and 5.4 of SMs.It is well-known that If |�P| ≪ |eR/�| , the uncertainty may not be a serious problem 82,83 .This condition, however, is not generally satisfied by all the compounds such as LiOsO 3 .Therefore, in Sect.5.4 of SMs, we have forced the transition to occur slowly from the starting structure = 0 " to the end structure = 1 " through the intermediate structures = 0.1, 0.2, ..., 0.9 .To this end, we have constructed the first intermediate structure = 0.1 " to be very close to the starting structure " = 0 ", as discussed in detail in Sect.5.4 of SMs.By comparing structures " = 0 " and " = 0.1 ", we have introduced some atomic vector steps 0.1 and distorted the structures one by one to gradually and slowly arrive at the endpoint " = 1 " step by step, see Sect.5.4 of SMs.In this way, we find a chance to identify a sudden change (jump), if any, in the calculated polarization at an intermediate distorted structure compared to its previous and next structures.If a jump (ascent or descent) occurs, we modify it to make smooth the path by shifting the jumped polarization, i.e. pulling downward the ascent polarization or pushing upwards the descent polarization, using a negative or positive integer multiple of the quantum of polarization, as practically discussed below.In fact, by this way, we unwrap the polarizations (Berry phases) of the constructed structures step by step which are by default traditionally wrapped into the interval [−eR/2�, eR/2�] ≡ [−enc/2�, enc/2�] ( [−π , π] ).Unwrapping refers to adjusting the phases of a signal to allow for smooth transitions.When phase jumps between successive signals are greater than or equal to the difference of π , unwrapping the phase helps in achieving continuous signals.
The results show that the electronic components ϕ el and φ el and thereby entirely shift them so that they start from zero.These shifts by constant values do not change the results, because the spontaneous polarization as the final important physical quantity is obtained from the difference between the polarizations calculated at the starting and end structures, �P = P ( =1) − P ( =0) , so that any constant shifts are canceled out.The shifted φ (I * ) el and φ el are shown as φ (I * ) el and φ el in Fig. 4b.As the parameter increases from 0 to 0.8, the ionic component ϕ ion decreases from zero to near the lower limit of the border shown by a horizontal dashed line at ϕ = −0.5 , or equivalently at ϕ = 2π ϕ = 2π(−0.5)= −π , see Fig. 4a.Then, ϕ ion at " = 0.9 " suddenly jumps to near the upper border indicated by a horizontal dashed line at ϕ = 0.5 , or equivalently at ϕ = 2π ϕ = 2π(0.5)= π , see Fig. 4a.Like before " = 0.9 ", again ϕ ion continues to decrease from " = 0.9 " to " = 1 ", see Fig. 4a.
The jump in ϕ ion detected at " = 0.8 " is not physically meaningful.By increasing the distortions very slowly, it may be expected to observe a smooth evolution between sequential structures leading to a non-zigzag path.Therefore, we unwrap, as shown in Fig. 4b, the sudden jump by pulling downwards the ascent ϕ ion at " = 0.9 " to ϕ ion which is performed by subtracting a quantum of Berry phase divided by 2π from ϕ ion , viz.
The ϕ (II) el remains unchanged, since ϕ (II) el needs to be neither unwrapped nor shifted.It is also represented as ϕ (II) el in Fig. 4b, keeping in mind that ϕ (II) el = ϕ (II) el .Although φ can be first unwrapped similar to ϕ ion → ϕ ion and then shifted similar to φ el → φ el , we obtain and represent it as φ in Fig. 4b more simply by the summation of the shifted φ el and wrapped ϕ ion as φ = ϕ ion φ el .The unwrapping and shifting procedure depicted in Fig. 4b yields smooth evaluations of the Berry phases across all structures-initial, intermediate, and final.
Analogous to the Berry phases shown in Fig. 4a, partial electronic and ionic components of the polarizations calculated in µC/cm 2 for the 11 structures " = 0, 0.1, 0.2, ..., 0.9, 1 " are shown in Fig. 4c.The polari- zations of the initial and final structures " = 0 " and " = 1 ", as the endpoints of the paths, are identical to the polarizations calculated in "Mean-field-like mBp method of polarization" and Sect.3.1 of SMs, compare Table 1 with the endpoints of the paths shown in Fig. 4c, taking the conversion relation "100 µC/cm 2 = 1 C/m 2 " into account.As expected from Eqs. ( 3) and ( 9) of SMs, the polarizations vary the same as Berry phases with respect to , compare Fig. 4a,c.Based on the wrapped results presented in Fig. 4c, the spontaneous polarization is calculated to be 94.55 µC/cm 2 , which is consistent with the results presented in Table 1 as 0.9455 C/m 2 , viz.�P = P ( =1) − P ( =0) = 94.55 µC/cm 2 = 0.9455 C/m 2 .In order uniquely determine P and make cer- tain the latter result, we determine the integer n by unwrapping the polarizations presented in Fig. 4c.To this end, the polarizations are unwrapped and the results are shown in Fig. 4d, where the tilde symbol in this figure denotes any necessary unwrapping and/or shifting tasks.The transformation procedure of the polarizations, including both unwrapping and shifting operations, from Fig. 4c to d is similar to that of the Berry phases from Fig. 4a to b, as discussed above.The variations of the unwrapped polarizations with respect to also behave like the unwrapped Berry phases, compare Fig. 4b,d and see the proportionality relations between polarizations and Berry phases in Eqs. ( 3) and ( 9) of SMs.The integer n is indicated by the vertical axes on the right of the Fig. 4c,d.By comparing Fig. 4c with Fig. 4d, the integer number n is determined to be −1 for the FE-LM LiOsO 3 .Therefore, by recalling that �P := P ( =1) − P ( =0) + enc/� and noting that ec/� = 71.78µC/cm 2 , ultimately we uniquely determined the spontaneous polarization of the FE-LM LiOsO 3 as � P = P ( =1) − P ( =0) = P ( =1) − P ( =0) + enc/� = 94.55 µC/cm 2 + (−1) × 71.78 µC/cm 2 = 22.77 µC/cm 2 .
The auxiliary symbols ∧ and ∼ are used only temporarily in this subsection for clarity.They indicate where wrapping and/or shifting operations are performed.It is important to note that wrapping and shifting are practical operations.Their sole purpose is to ascertain the final, factual spontaneous polarization.Spontaneous polarization is a physical quantity that can be observed in nature and measured experimentally.The theoretical calculation of spontaneous polarization may depend on these practical operations.However, the experimentally measured spontaneous polarization, as observed in nature, is evidently independent of these operations.Given this, we have chosen not to use the ∧ and ∼ symbols in other sections and subsections of this manuscript.We have simply reported the results without any additional symbols, such as P , P , and so on.This is under the understanding that the aforementioned operations are applied as necessary.

Validity of the SEP predicted for FE-LM LiOsO 3
In *"Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch", we predicted the SEP of ferroelectric lithium osmate (LiOsO 3 ) to be approximately 22.77 µC/cm 2 .This value aligns closely with the established SEP of Barium Titanate (BaTiO 3 ), reported to be 25 µC/cm 267,68 , as can be observed in Fig. 5a,c.This parallel is further corroborated by the research conducted by Zabalo et al. 30 , where similar critical bending radii were noted for these two compounds.Now, we aim to further validate and verify the accuracy of our prediction, aligning with predictions made by previous researchers as thoroughly reviewed in the introduction section.

Numerical verification: consistency of mBp and mWf
In "SEP of LiOsO 3 : mWf approach of electric polarization", the SEP of theFE-LM LiOsO 3 was numerically calculated to be 0.9496C/m 2 using the mWf method, as seen in Table 1.This value is known to suffer from the uncertainty problem, arising from the fact that a Wannier center is well-defined modulo R , as discussed in "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch".Analogously to the unwrapping procedure elaborated in detail in "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch" for the mBp method, we similarly found that n = −1 for the mWf method.This leads to � P = P ( =1) − P ( =0) = P ( =1) − P ( =0) + enc/� = 94.96µC/cm 2 + (−1) × 71.78 µC/cm 2 = 23.18µC/cm 2 for the SEP of the FE-LM LiOsO 3 by the mWf method, as shown in Fig. 4d and Table 3.By comparing the 23.18 µC/cm 2 predicted by the mWf method with the 22.77 µC/cm 2 predicted by the mBp method, we can at least conclude that the numerical mBp and mWf methods of calculating polarization consistently yield similar results for the SEP in this ferroelectric metal, as shown in Fig. 4d and Table 3.Furthermore, since the numerical mBp and mWf methods of calculating polarization, as discussed in "Mean-field-like mBp method of polarization", "mWf method of polarization" as well as Secs.3.1 and 3.2 of the SMs, are different, it is highly unlikely that their similar results were merely coincidental.However, in the following sections, we will provide further evidence to suggest that these results are also likely to be close to the experimental value.

Empirical verification: quadratic order
Here, we demonstrate that the spontaneous polarizations, as predicted by the mBp or mWf methods for the FE-LM LiOsO 3 , can be effectively fitted to the empirical equation proposed by Abrahams et al. 67 .The spontaneous polarization, P , and the phase transition temperature, T c , play vital roles in ferroelectrics 102 .These two funda- mental properties, i.e., P and T c , can be influenced by the atomic displacement (�z) of the homopolar atom, which is the most crucial quantity in "Displacive" ferroelectrics 102 .In compounds having similar symmetries, the homopolar atoms usually behave in a similar manner during ferroelectric phase transitions.Consequently, the variations in T c can be empirically estimated based on P 67,102,103 .
Here, we use empirical quadratic Eq. ( 12), derived directly from experimental data, and calculate the spontaneous polarization for the compound under investigation.To this end, we substitute the experimentally measured T c = 140 K for LiOsO 3 15 into Eq.(12).In this way, we empirically obtain �P = (21.50± 0.64) µC/cm 2 for the FE-LM LiOsO 3 (see Sect. 7.2 of SMs for details of the derivation).Our empirical result agrees with the theoretical predictions made using the mBp and mWf methods of polarization (see Table 3).This empirical evidence further strengthens the conclusion drawn at the end of "Numerical verification: consistency of mBp and mWf ".The consistency achieved based on quadratic Eq. ( 12), together with the agreement between results predicted by the mWf and mBp methods of polarization, helps to affirm the accuracy of our theoretical predictions.Not only does experimental data support Eq. ( 12) and, as a result, validate our numerical predictions, but the proportionality between T c and (�P) 2 also has strong theoretical backing, as discussed below.

Phenomenological verification: Landau-Ginzburg theory
Landau and Ginzburg developed a phenomenological theory for the second-order phase transition in ferroelectric materials by considering spontaneous polarization as an order parameter and expressing free energy with respect to this order parameter, P .The free energy density F P within this framework can be represented in the absence of external electric field and stress as follows 104,105 : where the power series is truncated at the fourth order, α 0 is a constant that depends on the materials, β is the Landau coefficient, and the critical phase transition temperature T c is the Curie-Weiss temperature T C .For the ( 13) Table 3. P calculated for LiOsO 3 , LiNbO 3 , LiTaO 3 and BiFeO 3 , together with the available experimental, empirical, and theoretical results.R is the magnitude of the real-space lattice vector along the polarization direction, i.e.R = c, where c is the lattice constant along the Cartesian z-axis.a is the other hexagonal lattice parameter, is the volume of the unit cell, eR is the quantum of polarization, and T c is the Curie temperature.Scheme, exchange-correlation functional (XC), Hubbard parameter (U), spin-polarization (SP), magnetic ordering (Order), Metallic state (Metal) and the Code are indicated.The units of the quantities are indicated.
The results calculated in the present work are denoted by *. a The datum is extracted from Ref. 67 , produced using an empirical equation proposed by Abrahams et al. b The lattice vector in [111] direction and the volume of rhombohedral unit cell, as extracted from Ref. 22 , are noteworthy for their distinct feature: the volume of the rhombohedral unit cell is one-third the size of a hexagonal structure.c The SEP predicted, in "Hypothetical verification using neural network: bandgap opening by imposing distortion", by the neural network at zero biaxial strain based on the data extracted from Ref. 27  we can deduce either P = 0 , representing the equilibrium polarization of the paraelectric phase, or ta nontrivial (nonzero) equilibrium polarization of the ferroelectric phase: Thus, the Landau-Ginzburg theory of second-order phase transition (paraelectric↔ferroelectric) implies, according to Eq. ( 15), that: We notice that the theoretical quadratic Eq. ( 16), which is derived from the Landau-Ginzburg theory of second order phase transition, is consistent with the empirical quadratic Eq. ( 12).As shown in "Empirical verification: quadratic order", our results are also consistent with the empirical Eq. ( 12), derived from experimental data.Consequently, our results are theoretically supported by Eq. ( 16) as well.

Empirical verification: biquadratic order
In "Phenomenological verification: Landau-Ginzburg theory", we truncated the power series of the free energy density F P (as expressed by Eq. ( 13)) at the fourth order.This resulted in a quadratic order polynomial for the transition temperature in terms of spontaneous polarization (see Eq. ( 16)).In "Empirical verification: quadratic order", we used Eq. ( 12), which also represents a quadratic equation.However, to incorporate a broader range of interactions, including anharmonic vibrations, it would be necessary to involve higher orders of spontaneous polarization.This approach can provide a more precise tool for verifying the reliability of our SEP calculations for FE-LM LiOsO 3 .With this in mind, let us go beyond the quadratic order.In principle, the power series of the free energy density could be truncated at the sixth order.This would, in principle, lead to a quartic (biquadratic) relation between spontaneous polarization and the transition temperature.Following this approach, we derive empirical biquadratic polynomials using available experimental P and T c data for a wide range of ferroelectric compounds.
In a T c − (�P) 2 Cartesian coordinate system, as shown in Fig. 5a, we have positioned thirteen points (T c , (�P) 2 ) .Twelve of these points are based on experimental data 68,85,87,89,92-101 for various FE perovskite semi- conductors and one of them, whose experimental T c is taken from Ref. 15 , is our numerical data calculated by the mBp method for the FE metal under question.In this figure, we did not include the numerical point calculated by the mWf method because it closely aligns with that calculated by the mBp method for the FE LiOsO 3 metal.We only fitted a quartic-order polynomial to ten of the thirteen experimental points.These 10 experimental points are indicated by hollow-square-symbols, see Fig. 5a.Some of the experimental data 85,89,97 have been recently measured.Fig. 5a includes two different experimental data points for LiTaO 3 .We used one of these data points, recently measured by Zhang et al. 89 , for a better fit, as it brings the coefficient of determination (R-squared) closer to unity, see R 2 = 0.97 in Fig. 5a.Similarly, for NKbO 3 compound, among two different sets of experimental data included in Fig. 5a, we used only one of them, reported by Günter 96 , for a more efficient fit that brings R 2 closer to unity.We have excluded our numerical point, indicated by solid-square-symbol, in the fitting procedure to avoid affecting the resultant empirical quartic relation emerged from the experimental data.By this way, we have obtained the following biquadratic empirical relation: where a linear behavior of T c versus (�P) 2 , in agreement with Eq. ( 12) reported by Abrahams et al. 67 , can be observed for smaller P for which the effects of the quartic term due to the factor of 10 −5 can be less than the quadratic term, see Fig. 5a.However, a deviation from the linear behavior of T c with respect to (�P) 2 can be observed for larger P , where the quartic term can compete against the quadratic term despite the 10 −5 factor.Even though the curve shown in Fig. 5a is fitted to the experimental data, and the point (T c , (�P) 2 ) correspond- ing to LiOsO 3 is absent, the fitted curve closely aligns with our numerical results calculated for LiOsO 3 .This validates our spontaneous polarization calculated for the FE-LM LiOsO 3 .
Additionally, we have limited the number of compounds utilized in the fitting procedure to the experimental data in Fig. 5a to just LiTaO 3 and LiNbO 3 , which belong to the same family as LiOsO 3 displayed in Fig. 5b.These LiXO 3 (X=Nb, Ta) ferroelectrics, similar to LiOsO 3 , undergo an R3c to R 3 c phase transition and possess symmetry comparable to FE-LM LiOsO 3 .Within these three compounds-LiNaO 3 , LiTaO 3 , and LiOsO 3 -their homopolar atoms (Ta, Nb, and Os) play key roles in determining their ferroelectric properties.( 14) www.nature.com/scientificreports/A quadratic polynomial was fitted to the experimental (T c , (�P) 2 ) points considering only LiTaO 3 and LiNbO 3 , resulting in the following biquadratic polynomial for T c as a function of P: Its coefficient of determination, R 2 = 0.99 , is closer to unity than R 2 = 0.97 , determined for Eq.( 17), indicating a better fit to the experimental data.To further test the accuracy of our results, LiOsO 3 was also considered.The data for LiOsO 3 were not used to obtain the fitted curve in Fig. 5b.However, this curve still closely approximates the point associated with LiOsO 3 , lending credibility to our theoretical predictions.
Lebeuglea et al. reported a large SEP for the multiferroic BiFeO 3 compound 85 .In Fig. 5c, we have included this multiferroic compound along with the standard ferroelectric compounds from Fig. 5a.The data from Lebeuglea et al. 85 for BiFeO 3 were added to Fig. 5c, and a polynomial was fitted to the expanded dataset.This process yielded the following biquadratic polynomial with R 2 = 0.98 for T c as a function of P: This empirical polynomial also validates our calculated polarization for the LiOsO 3 , see how the point (T c , (�P) 2 ) of LiOsO 3 is close to the curve in Fig. 5c.

Hypothetical verification: bandgap opening by Hubbard model using GGA+U
In this subsection, our objective is to calculate the SEP of LiOsO 3 in its hypothetical semiconductor phase.To do this, we slightly open the bandgap using PBE-GGA+U.This adjustment enables us to use the conventional Bp or Wf method without any restrictions or modifications.Our goal is to demonstrate that the SEP, as projected by the conventional Bp or Wf method for the synthetically developed G-AFM semiconductor LiOsO 3 , is comparable to the SEPs predicted by the mBp or mWf polarization method for the naturally occurring nonmagnetic metal LiOsO 3 .This comparison is intended to validate that our modified mBp or mWf method aligns well with the conventional Bp method, even without our modifications.
To achieve our aim, we seek to validate that the drawn conclusions are robust against changes in DFT methods and supercell configurations.Hence, we will proceed with the calculation of SEP in a manner that harmonizes more consistently with the methodologies used in this research.Thus, following our study's approach, we will also estimate the SEP using the conventional Berry phase theory, as implemented in the WIEN2k code based on full potential.To this end, we will also utilize GGA+U to open the gap and transform the rhombohedral unit cell into a hexagonal supercell of LiOsO 3 , as detailed in Secs.5.2 and 5.3 of the SMs, while imposing G-type AFM ordering.In our GGA+U calculations, we will set U eff = 0.2 eV .This choice allows us to open the bandgaps of both polar noncentrosymmetric R3c and non-polar centrosymmetric R 3 c structures, by minimal amounts of 0.061 eV and 0.071 eV respectively.Indeed, the aim of opening the gap was to leverage the standard Berry phase of polarization to corroborate our findings.By applying GGA+U with U eff = 0.2 eV , we computed the electronic, ionic, and total Berry phases for the initial, final, and the nine intermediate superstructures, as elaborated in Sect.5.4 and illustrated in Fig. SM6 of the SMs.
As illustrated in Fig. 6a, the electronic Berry phase exhibits a linear trend, while both ionic and total Berry phases present a zigzag pattern across the array of 11 superstructures.The abrupt zigzag fluctuations in the ionic Berry phases are tempered by unwrapping sudden shifts, thus creating a smoother line, as demonstrated in Fig. 6b.We also modify the trajectory of the electronic Berry phase to ensure its origin at " = 0 " aligns with the zero-point, as shown in Fig. 6b.The total Berry phase is then obtained by combining the electronic and ionic Berry phases (Fig. 6b).
The initial or misplaced polarizations, visible in Fig. 6c, are adjusted as depicted in Fig. 6d.This specific method and its aim are extensively detailed in "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch" and visually illustrated in Fig. 4, so further elaboration here is redundant.In Fig. 6, we have intentionally omitted auxiliary symbols ∧ and ∼ , as discussed in "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch".Upon comparing Fig. 6c and d, we find that the SEP, calculated by considering only the initial " = 0 " and the final " = 1 " structures, aligns with the SEP computed across all 11 structures, evidenced by the identical value of �P = 24.33 µC/cm 2 .Hence, the consideration of intermediate superstructures, or adjustment of abrupt zigzag patterns, is redundant for this compound.However, this finding should not be generalized, as the SEP may vary when intermediate structures are included, as demonstrated in Fig. 4c,d.
Interestingly, the phases (polarizations) for the G-AFM phase fall within the range [−2π, 2π Rad]([−71.78,71.78 µC/cm 2 ]) (Fig. 6a,b), while the phases (polarizations) for the NM phase are confined to [−π , πRad]([−71.78/2, 71.78/2 µC/cm 2 ]) (Fig. 4a,b).This discrepancy stems from spin-polarized calculations for the G-AFM phase (Fig. 6a,b) versus non-spin-polarized calculations for the NM phase (Fig. 4a,b).In spin-polarized calculations, the Berry phases (and corresponding polarizations) for spin-up states are summed with those for spin-down states, causing the doubling in wrapping intervals due to the additional spin degree of freedom.However, this factor does not affect the calculation of P , where the difference between initial and final polarizations is calculated.Our computed SEP, �P = 24.33 µC/cm 2 for the nonmetallic state, aligns with the SEP calculated for the metallic state of LiOsO 3 , �P = 22.77 µC/cm 2 , as seen in Fig. 4d and Table 3, along with Fig. 6d.
Continuing our analysis, we acknowledge the work of He et al., who calculated the SEP as 22.23 µC/cm 222 .Their research relied on the standard Berry phase theory 52,55 and employed the pseudopotential-based VASP code 46-50 .Using LSDA+U with U eff = 2 eV , they established a 2 × 2 × 2 supercell derived from the rhombohedral unit cell of G-AFM LiOsO 3 .
The computed SEP of 22.23 µC/cm 2 for the conjectured G-AFM insulating LiOsO 3 aligns with our SEP for the authentic nonmagnetic metal LiOsO3, calculated as 22.27 (23.18) µC/cm 2 by the mBp (mWf) method of polarization we have developed (see Fig. 4d and Table 3).This compatibility validates our mBp and mWf methods and indicates that the influence of GGA+U on the class I * , where the valence bands cross the Fermi level, is limited (refer to Fig. 1 and Table 1, where el is reported as a moderate value of 0.0746 C/cm 2 for the polar NCS R3c structure).
Before concluding this subsection, we will examine the implications of U on polarization in greater depth.Additionally, see Sect.5.6 of the SMs for a discussion regarding the influence of U on lattice parameters and the energy bandgap.
Our investigation reveals a marginal discrepancy between the SEP obtained using GGA+U with U eff = 0.2 eV and that determined by He et al. 22 using LDA+U with U eff = 2.0 eV , where �P = 24.33 > 22.23 µC/cm 2 (refer to Table 3).Interestingly, the SEP value acquired by He et al. aligns more closely with our PBE-GGA calculated outcome for the metallic state than with our GGA+U calculations involving U eff = 0.2 eV , such that �P = 22.77 ≈ 22.23 µC/cm.
In alignment with the procedure employed by He et al., we recalculated the SEP using GGA+U with U eff = 2.0 eV .The resultant SEP value, �P = 22.32 µC/cm , displays superior concurrence with the findings reported by He et al., as well as with our mBp and mWf polarization-derived SEPs in the metallic state.These shared findings indicate �P = 22.32 ≈ 22.77 ≈ 22.23 µC/cm (refer to Table 3 and Fig. 4d).
Nonetheless, it is essential to remember that we applied the PBE-GGA functional and the FP-APW+lo method, differing from those used in Ref. 22 .Thus, although the closer agreement when using the identical U value is not surprising, given these divergences in computational methodology, it warrants mention.
Given the above analysis, the effects of U parameter on the SEP are significant, albeit not necessarily large.Therefore, we have systematically studied these effects, as promised earlier.We have calculated P by GGA+U for a range of U from 0.2 to 2.4 eV in steps of 0.2 eV, as shown in Fig. 7a.Our findings reveal a decrease in P with increasing U.
This reduction is attributable to the electronic part of the polarization of the polar NCS R3c phase, which corresponds to the final structure " = 1 ".For clarification, let's examine the changes in the electronic and ionic parts of the Berry phases in both the initial " = 0 " and final " = 1 " structures.
The total polarizations for the initial structure " = 0 ", represented by blue empty square symbols in Fig. 7b, show fluctuations between −ec/2� and +ec/2� .By subtracting a quantum of polarization, +ec/� , from P ( =0) = ec/2� at selected U values, we smooth the path of polarization.
Therefore, P ( =0) = −ec/2� remains constant for U ∈ [0.2, 2.4 eV] , as indicated by orange empty circle symbols in Fig. 7b.This indicates that P ( =0) remains constant at half of the quantum of polarization and is not affected by variations in U for U ∈ [0.2, 2.4eV] .Hence, the variation of P cannot be attributed to the non-polar CS R 3 c structure and originates solely from the polar NCS R3c structure.
Our GGA+U results indicate that P ( =1) ion remains constant at −0.40482 µC/cm 2 for U ∈ [0.2, 2.4 eV] , while P ( =1) el decreases as U increases, as depicted in Fig. 7d.Consequently, the variation of P with U is solely attribut- able to the variation of P ( =1) el with U, as shown by the comparison of Fig. 7a with Fig. 7c.The band structures calculated by PBE-GGA for the NM phase and by GGA+U with U = 0.2 eV and U = 2 eV for the G-AFM phase, as presented in Fig. SM1 of the SMs, reveal 18 bands within the energy range of [−-2, 2 eV].Our calculations show minimal changes to the valence bands with increasing U in GGA+U, while the conduction bands move away from the Fermi level.Our findings suggest that the variation in SEP is primarily due to changes in the electronic part of the NCS phase.
In conclusion, our study articulates the validity of SEPs as projected by the mBp and mWf methodologies for the NM metallic phase of LiOsO 3 .We have demonstrated that these SEPs can be accurately approximated and corroborated using the generalized gradient approximation plus Hubbard U (GGA+U) for the contrived G-AFM nonmetallic phase of LiOsO 3 .Furthermore, it has been discerned that the influence of diminutive U parameters on the SEP is relatively marginal, permitting us to open the bandgap of the G-AFM phase through GGA+U and estimate the magnitude of SEP for the NM metallic LiOsO 3 .Therefore, our endeavors have successfully met the stated objectives.

Hypothetical verification using neural network: bandgap opening by imposing distortion
Here, we present an alternative strategy for bandgap modulation, complementing the GGA+U methodology described in "Hypothetical verification: bandgap opening by Hubbard model using GGA+U".It involves the application of appropriate strain to the system.Zhang et al. 27 examined the correlation between the bandgap and biaxial strain in LiOsO 3 , maintaining consistent crystal symmetry.Their results demonstrated a direct relationship between the increase in the percentage of tensile biaxial strain, represented as ε (%) , and the augmentation of the bandgap.
However, the limitations of linear extrapolation became apparent when comparing the rate of change of P at different ε values.A quadratic fit was applied, resulting in the equation �P = [(−1.88± 0.57)ε 2 + (9.57± 4.55)ε + (2.10 ± 8.91)]µC/cm 2 , as seen in Fig. 8b.Despite an improved R 2 value for the quadratic function, the relative error remained significantly high.
Given these limitations, the study turned to machine learning, employing the Bayesian regularization-trained multilayer perceptron (MLP) methodology [106][107][108][109][110][111] .The MLP procedure applied Bayesian activation function to optimize the weights and reduce the error function E(O) = T − f (IW io ) .The network output vector O is computed as O = f (IW io ) .After the training, the predicted SEP at ε = 0 was 23.46 µC/cm 2 , as shown in Fig. 8c, consistent with the predictions made by the mBp and mWf methods.Lastly, a refined empirical quadratic equation was derived for P as a function of strain ε by including the prediction from the neural network, yielding: which was graphed in Fig. 8d.This model's prediction at ε = 0 aligns closely with the neural network prediction, as well as the predictions made by the mBp and mWf methods, providing validation for these methodologies.

Comparative verification: BiFeO 3 , LiNbO 3 , and LiTaO 3
In this section, we aim to validate the efficacy of our unwrapping procedure in determining the most suitable branch for the unique prediction of the resultant SEP, as elucidated and applied to LiOsO 3 in *"Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch" for its actual metallic state, and in "Hypothetical verification: bandgap opening by Hubbard model using GGA+U" and "Hypothetical verification using neural network: bandgap opening by imposing distortion" for its conjectural nonmetallic state.
To facilitate this, we have selected the normal ferroelectrics BiFeO 3 , LiNbO 3 , and LiTaO 3 for comparison.These compounds, similar to LiOsO 3 , undergo the R 3c-R3c ferroelectric transition.This transition, as detailed in Sect.5.4 of the SMs, facilitates the creation of their intermediate structures in a manner analogous to LiOsO 3 .We simulate this transition from = 0 to = 1 , passing through the intermediate structures = 0.1, 0.2, ..., 0.9 as illustrated in Fig. SM6 of the SMs.
Contrary to LiOsO 3 , the SEPs of these selected compounds have been both experimentally measured 85,87,89 and theoretically computed 84,88 .The availability of experimental data provides us the opportunity to evaluate the accuracy of our unwrapping and uniquifying process in the context of real-world results.
We begin our analysis with the multiferroic BiFeO 3 , drawing extensively from our previous work 78 .For its G-AFM phase, we calculate the Berry phases and polarizations across all 11 superstructures depicted in Fig. SM6 of the SMs.These calculations, conducted via GGA+U methodology with an optimized U eff = 4 eV , are presented in Fig. 9. Figure 9a presents the total and wrapped partial electronic and ionic components of the Berry phases, scaled by 2π.
The scaled wrapping interval for G-AFM BiFeO 3 , [−0.5, 0.5] , is precisely half of that for G-AFM LiOsO 3 , [−1, 1] , as observed when comparing Fig. 9a with Fig. 6a.Spin-orbit coupling (SOC) demonstrates a greater .The unit of polarizations is µC/cm 2 in (c) and (d).Like Fig. 6 but unlike Fig. 4, here, only for simplicity, the auxiliary symbols ∧ and ∼ are not used.All the Berry phases and as a result polarizations are calculated by the standard Berry phase scheme including SP and SOC by GGA+U with U eff = 4 eV for the Multiferroic G-AFM BiFeO 3 along the distortion path as functions of structure from " = 0 " to " = 1 " by step 0. influence on G-AFM BiFeO 3 than G-AFM LiOsO 3 , thus our computations for G-AFM BiFeO 3 (G-AFM LiOsO 3 ) incorporate (exclude) SOC, as illustrated in Fig. 9 (Fig. 6).
Relativistic quantum mechanics elucidates that the DOSs of spins up and down either combine or separate in accordance with the presence or absence of SOC [112][113][114] .Consequently, Berry phases, polarizations, and other electronic properties either couple or decouple accordingly, necessitating the choice of a relativistic ( |j, m j , l, s� ) or non-relativistic basis ( |l, m l , s, m s �).
In this context, the wrapping interval of Fig. 6a is twice that of Fig. 9a, given the corresponding quanta of polarizations.Fig. 9c illustrates the changes in the wrapped polarizations across the 11 superstructures.
The uncertainty problem arising from the calculated wrapped SEP of 16.98 µC/cm 2 , which differs significantly from the experimental value of 100.00 µC/cm 284 , requires resolution via identification of the optimal branch, or integer quantum number n.
We apply an unwrapping procedure (analogous to that discussed in "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch") to the Berry phases and polarizations, as depicted in Fig. 9b,d, respectively.By comparing total polarizations at " = 1 " before and after wrapping, or equivalently total Berry phases, we discern the integer quantum number n to be −2.
Therefore, taking into account n = −2 and the quantum of polarization 59.44 µC/cm 2 , the SEP derived from the unwrapped total polarization 16.98 µC/cm 2 is calculated to be −101.90µC/cm 2 .
Our SEP calculation using the density functional and standard Berry phase theories with GGA+U and U eff = 4 eV is consistent with both experimental data 85 and theoretical results 84 .This congruence affirms the validity of our unwrapping and uniquifying procedure.We also replicated standard Berry phase calculations using GGA+U with U eff = 2 eV , and found an increase in | P| as U decreases from 4 to 2 eV, affirming the P -U relationship depicted in Fig. 7a.
In this investigation, we consider the typical ferroelectrics (FEs) LiNbO 3 and LiTaO 3 .These materials exhibit equivalent electronic configurations, similar chemical compositions, and undergo identical R 3c-R3c ferroelectric transitions as LiOsO 3 84,86 .Their spontaneous electrical polarizations (SEPs) have been documented experimentally.In Fig. 10a-d, both the wrapped and unwrapped or shifted total and partial polarizations as a function of superstructures = 0, 1, 2, . . ., 1 are represented for LiNbO 3 and LiTaO 3 , respectively.The unit of polarizations is µC/cm 2 .Like Figs. 6 and 9 but unlike Fig. 4, here, only for simplicity, the auxiliary symbols ∧ and ∼ are not used.All the Berry phases and as a result polarizations are calculated by the standard Berry phase scheme including SP by PBE-GGA for the normal FE LiNbO 3 and by LDA for the normal FE LiTaO 3 along the distortion path as functions of structure from " = 0 " to " = 1 " by step 0.1.For comparison, experimental 87,87,89 and theoretical 84,86,88 SEPs are included in (b) for LiNbO 3 (in (d) for LiTaO 3 ).
These manipulations do not modify the SEPs for LiNbO 3 and LiTaO 3 ; for both, the shift n is null, demonstrated by an equivalent 78.31 µC/cm 2 value for LiNbO 3 and 53.12 µC/cm 2 for LiTaO 3 , observed in Fig. 10a-d for P.
The G-AFM FE LiOsO 3 , with a bandgap, also exhibits n = 0 , illustrated by a constant 24.33 µC/cm 2 value for P in Fig. 6c,d 84,86 .Conversely, the NM FE-LM LiOsO 3 and the GAM normal multiferroic BiFeO 3 dis- play n = 1 � = 0 and n = −2 � = 0 , respectively.For these cases, the SEPs differ as shown by the distinct values extracted from the respective figures for P .This indicates that the unwrapping or shifting procedure modifies " P ( =1) − P ( =0) " unless the optimal branch is chosen, in which case it remains consistent.
For LiNbO 3 (Fig. 10a), some discontinuities are noted, while for LiTaO 3 (Fig. 10c), these are absent.This suggests a n = 0 for LiNbO 3 , as the difference " P ( =1) − P ( =0) " is impervious to a simple shift.However, this is not solely dependent on the initial and final structures " = 0 " and " = 1 ", and necessitates examination of intermediate structures to guarantee the selection of the optimal branch.Theoretical data and experimental values at room temperature for LiNbO 3 are presented in Fig. 10b and Table 3, offering a comparative analysis.Our results align with existing theoretical research and experimental data 84,86,87 .Likewise, for LiTaO 3 , our findings accord with the experimental SEPs measured by Wemple et al. 87 , and Zhang et al. 89 , and the theoretical datum calculated by Tan et al. 88 .Notably, these consistencies between theoretical and experimental results extend further when considering that the theoretical DFT results were computed at zero temperature while experimental data were obtained at room temperature.Reduction in temperature can decrease entropy and enhance the electrical order, thus potentially increasing the experimentally measured polarization.This, in turn, can improve the consistency between theoretical results (calculated at zero temperature) and experimental data (recorded at lower temperatures).
The ferroelectric structural distortion in LiXO 3 perovskites ( X = Nb, Ta, Os ) has been linked to the hybridiza- tion of O:p and X:d orbitals 22,84,115 .The degree of this distortion is influenced by the eccentricity of the X atom, as measured by the c/a ratio, a recognized indicator of structural distortion strength 84 , and outlined in "SEP direction in LiOsO 3 FE-LM".
In the LiXO 3 family, the computed c/a ratios are 2.69, 2.67, and 2.64 for LiNbO 3 , LiTaO 3 , and LiOsO 3 respectively, illustrating a decrease in structural distortion with decreasing c/a ratio.This trend also corresponds to a decrease in the spontaneous electric polarization (SEP).Hence, the smaller SEP of LiOsO 3 can be substantiated when compared to the larger SEPs of LiNbO  3.
However, this comparison should not be generalized for compounds with different correlations.While the SEP decreases with an increase in U for LiOsO 3 (with a constant c/a ratio, as depicted in Fig. 7), it may increase with U if the c/a ratio is not fixed due to its own increase with U (as shown in Table (SM2) of the SMs).Hence, the SEP is influenced by both the c/a ratio and U, and its behavior relative to U is contingent upon whether the c/a ratio is held constant.This necessitates caution when extrapolating these results to other cases.

Conclusion
In this work, we have explored a groundbreaking possibility: the existence of nonzero spontaneous electric polarization (SEP) in metals.We have challenged the widely accepted belief that itinerant electrons invariably annihilate ferroelectricity in metals.Our work builds on the theoretical conjecture of Anderson and Blount, and experimental findings of Y. Shi and team, offering a quantitative validation.Our research's pivot is the adjustment of existing methods, namely the Berry phase theory and Wannier functions theory, for calculating electric polarization in systems with nonzero bandgaps.By addressing their limitations and modifying these methods, we have developed the modified Berry phase (mBp) theory of polarization and modified Wannier functions (mWf) theory of polarization.These adapted methods are poised to work effectively in predicting the SEP of metals.
In the case study of the ferroelectric-like metal (FE-LM), lithium osmate (LiOsO 3 ), our calculated SEP demonstrates an alignment with the SEP in Barium Titanate (BaTiO 3 ), a regular ferroelectric compound.The consistency with both empirical and theoretical data underscores the validity of the mBp and mWf methods.We further validated our findings via multiple approaches: numerical verification, empirical testing, comparison with the Landau-Ginzburg theory, hypothetical adjustments to the bandgap, and comparative analysis with normal ferroelectric materials.
Notably, due to the absence of SEP at zero biaxial strain, we employed the multilayer perceptron method -a subset of the feedforward artificial neural network class in machine learning -to project the SEP at this state.This prediction was based on available data for nonzero strains.
Overall, the consistency across all validation methods implies a high likelihood of our predictions being accurate.As such, our proposed mBp and mWf methods can be reliably applied to predict SEP in metals.This opens up avenues for the theoretical identification and practical synthesis of new ferroelectric-like metals, enhancing the broader understanding of ferroelectricity, as aptly expressed by Evgeny Y. Tsymbal.Our research marks a significant contribution to the continuous enrichment of knowledge in physics and material science, 100 years after the original discovery of ferroelectricity.

Figure 2 .
Figure 2. The first Brillouin zone used for the calculations of the electronic part of the Berry phase in the reciprocal lattice, decomposing k-points into k and k ⊥ samples.Here, k is parallel to and k ⊥ is perpendicular on the polarization direction G µ , as shown on the right side of the figure.Several 2D-planes are formed by the k ⊥ with the normal vectors k .The points k are distributed over the parallel strings.Occupation numbers of the bands at k are indicated by n(k ) .The contribution of the electronic Berry phase for the structure at k are shown by ϕ ( ) el,µ (k ) .By the formula, as indicated in the top-left side of the figure, and considering ϕ ( ) el,µ (k ) , included in the overlap integral O ( ) M×M (k j , k j+1 ) , the Berry phases are calculated along each k strings individually, and the results for the structure are mapped as ϕ ( ) el,µ (k ⊥ ) into the central 2D-sheet with the area A ⊥ .Occupation numbers of the bands at k ⊥ are indicated by n(k ⊥ ) By the formula, as indicated in the bottom- left side of the figure, the average of the mapped Berry phases are calculated over the area A ⊥ .For more detail see Sect.1.1 of SMs.
(I * ) el ) is the partial electronic contribution of Wannier centers of class II (I * ) into the electronic part of polarization P ( ) which gives the occupancy of the corresponding center of charge.The values of these areas are the desired occupation numbers n W ( ) n,R .In the fourth step, we multiply each of the Wannier centers r W ( ) n,R obtained in the second step by their corresponding occupation numbers n W ( ) n,R obtained in the third step.By this way, n W ( )

(
),(I * ) el are calculated.Now, by adding P ( ),(I * ) el calculated in this step to P ( ),(II) el calculated in the first step, we obtain the electronic polarization P ( )

Figure 3 .
Figure 3. (a) DOS projected on one of the 6 Wannier centers related to the energy interval I, as shown in Fig. 1a.The initial DFT calculations were performed by the PBE-GGA DFT for the LiOsO 3 in its polar NCS rhombohedral structure.The colored area under the projected DOS up to the Fermi level shows the occupancy of the Wannier center n at R for the polar NCS rhombohedral structure of the compound, i.e. n W ( ) n,R .The area is calculated by taking integration of the projected DOS up to the Fermi level.(b) The real-space plot of maximally localized functions, W ( ) n,R , constructed from the Bloch states.(c) Original bands were generated directly from the PBE-GGA DFT calculations for the polar NCS phase of the rhombohedral LiOsO 3 , see thin black bands.Wannier-interpolated bands obtained from the subspace selected by an initially unconstrained projection onto atomic Os:d z 2 , Os:d x 2 −y 2 , Os:d xy orbitals for the isolated class of bands I, see thick blue bands.The Fermi level is set to zero in both the DOS and band structure figures.
(II) el align with the c-axes of the hexagonal supercells in both the CS R 3 c and NCS R3c phases.Besides the polarization directions, the computed values of P ion for both the centrosymmetric (CS) R 3 c and non-centrosymmetric (NCS) R3c phases, as well as P (II) 3.2 of the SMs, we focus on energy interval I to calculate the positions of the Wannier centers.These positions are determined by evaluating the integral r|W ( ) n,R (r)| 2 dr , which results in the position vector Vol.:(0123456789)Scientific Reports | (2024) 14:672 | https://doi.org/10.1038/s41598-023-49463-wwww.nature.com/scientificreports/�r�Wn, R ( ) = (�x�Wn, R ( ) , �y�Wn, R ( ) , �z�Wn, R ( ) ) .These calculations are performed for each Wannier center, numbered from n = 1 to J (I) , while considering both the CS R 3 c and NCS R3c phases individually.The number of bands of class I, M (I) , is 18 which equals the number of Wannier centers of class I, i.e.J (I) = 18 , for the hexagonal supercells of both the CS R 3 c and NCS R3c phases.In the FP-LAPW DFT calculations, we have set the separation energy of the valence electrons from core electrons to −9.0 Ry , leading to 306 valence electrons.The number of fully occupied bands of class II is 144.These bands contain 288 = 144 × 2 electrons.Thus, the bands of class I contain 18 = 306 − 288 valence electrons.These 18 electrons mostly come form the Os atoms, see Fig. SM2 of SMs.Although the DOSs shown in Fig.
SM2 of SMs.By considering these 3 valence d-electrons, it can be also verified that the total number of valence electrons of class I are approximately 3 × 6 = 18 .These 18 valence electrons are distributed over the bands of class I, including valence and low-lying conduction bands, so that 6 bands (containing ≈ 12 electrons) is almost fully occupied, 6 bands (containing ≈ 6 electrons) is partially occupied, and 6 bands (containing ≈ 0 electrons) remain almost empty.In fact, these 18 valence electrons are distributed over the bands of class I * , including valence bands only.The polarizations P ( ),(I * ) el

2 AAB 3 AAB 4 AAB 5 AAB 6
leads to18.00 for both the CS R 3 c and NCS R3c phases individually.It is worth noting that the number of Wannier cent- ers of class I, J (I) , in J (I) =18 n=1 n W n,R = 18.00 equals the number of bands of class I, M (I) , while the resultant value of the summation yields 18.00 which equals the number of valence electrons of class I * .This verifies that the occupation numbers are correctly calculated and the mWf procedure works well so far up to this step.The results show that the occupation numbers of the Wannier centers n W ( ) n,R are almost either 1.09 (:= A) or 0.82 (:= B) which are close to unity, viz.A = 1.09 ≈ B = 0.82 ≈ 1.00 .This shows that the occupation numbers can be approximately halved by including low-lying empty conduction states besides the fully occupied valence states for constructing the maximally localized Wannier centers.This shows that the 18 electrons are almost uniformly distributed over the 18 centers of the maximally localized Wannier functions constructed from both valence and conduction states.More precisely, by taking the differences between the values of A = 1.09 and B = 0.82 into account, a sequence AAB 1 AAB with a repeating pattern AAB involving 3 elements can be observed which is periodically repeated 6 times for both of the phases.If we multiply the number of elements of the repeating pattern, 3, by the number of repetitions of the pattern, 6, we obtain the number of the 18 centers associated with the maximally localized Wannier functions constructed from both valence and conduction states, viz. 3 × 6 = 18 .Approximately half of these states, ≈ 9 , belong to the valence region and the other half belong

Figure 4 .
Figure 4. (a) Total, and wrapped partial Berry phases versus .(b) Unwrapped total, and partial Berry phases versus .All the Berry phases are scaled by 2π in (a,b) so that the interval of wrapping is simplified from [−π , π] to [− 0.5, 0.5] in (a).(c) Total, wrapped partial, and corresponding spontaneous polarizations versus .The partial polarizations are wrapped into [ −ec/2�, ec/2� ], where ec/� = 71.78µC/cm 2 is the quantum of polarization.(d) Unwrapped total, partial, and corresponding spontaneous polarizations versus .The unit of polarizations is µC/cm 2 in (c,d).The auxiliary symbols ∧ and ∼ indicate that when wrapping and/or shifting are/is performed, if necessary, compared to the results presented in Table1, see "Uniquifying of spontaneous polarization of LiOsO 3 by finding the best branch" where the symbols are defined.All the Berry phases and as a result polarizations are calculated by the mBp scheme including non-spin-polarized PBE-GGA along the distortion path as functions of structure from " = 0 " to " = 1 " by step 0.1.The quantum of polarization and its number n are shown in (c,d).Our SEPs calculated by mWf, PBE-GGA+U with U = 0.2 and 2 eV are presented for comparison.The mBp, mWf, and empirical results presented in this figure are obtained for the metallic state of the NM LiOsO 3 .Our GGA+U results and the GGA+U result taken from Ref.22 are calculated for the nonmetallic state of the G-AFM LiOsO 3 using the standard Berry phase method.The P Emp is extracted from Ref.67 , generated using an empirical equation.The theoretical datum is taken from Ref.22 .

Figure 6 .
Figure 6.(a) Total, and wrapped partial Berry phases versus .(b) Unwrapped total, and partial Berry phases versus .All the Berry phases are divided by 2π in (a) and (b) so that the wrapping interval is converted from [−2π, 2π] to [− 1, 1] in (a).(c) Total, wrapped partial, and corresponding spontaneous polarizations versus .The partial polarizations are wrapped into [ −ec/�, ec/� ],where ec/� = 71.78µC/cm 2 is the quantum of polarization.(d) Unwrapped total, partial, and corresponding spontaneous polarizations versus .The unit of polarizations is µC/cm 2 in (c) and (d).Unlike Fig.4, here, only for simplicity, the auxiliary symbols ∧ and ∼ are not used.All the Berry phases and as a result polarizations are calculated by the standard Berry phase scheme including SP using GGA+U with U = 0.2 eV for the nonmetallic state of the G-AFM LiOsO 3 along the distortion path as functions of structure from " = 0 " to " = 1 " by step 0.1.The quantum of polarization and its number n are shown in (c) and (d).

Figure 7 .
Figure 7. (a) Spontaneous polarization, P , (b) total polarization of the CS phase, P ( =0) , (c) total polarization of the NCS phase, P ( =1) , (d) electronic polarization of the NCS phase, P ( =1) el , and (e) ionic polarization of the NCS phase, P ( =1) ion , calculated by the standard Berry phase approach for the nonmetallic state of the G-AFM LiOsO 3 using GGA+U with U ∈ [0.2, 2.4 eV] .The jumps in P ( =0) occurred at U = 0.6, 1.0, 1.4, 1.6, 1.8, 2.0 , as shown by empty blue square symbols in the inset (b), are pulled down and the blue zigzag path is made smooth as a straight orange line by subtracting one quantum of polarization, see empty orange circle symbols in the inset (b).

( 21 )Figure 9 .
Figure 9. (a) Total, and wrapped partial Berry phases versus .(b) Unwrapped total, and partial Berry phases versus .All the Berry phases are divided by 2π in (a) and (b) so that the wrapping interval is converted from [−π , π] to [− 0.5, 0.5] in (a).(c) Total, wrapped partial, and corresponding spontaneous polarizations versus .The partial polarizations are wrapped into [ −ec/2�, ec/2� ], where ec/� = 59.44 µC/cm 2 is the quantum of polarization.(d) Unwrapped total, partial, and corresponding spontaneous polarizations versus.The unit of polarizations is µC/cm 2 in (c) and (d).Like Fig.6but unlike Fig.4, here, only for simplicity, the auxiliary symbols ∧ and ∼ are not used.All the Berry phases and as a result polarizations are calculated by the standard Berry phase scheme including SP and SOC by GGA+U with U eff = 4 eV for the Multiferroic G-AFM BiFeO 3 along the distortion path as functions of structure from " = 0 " to " = 1 " by step 0.1.The quantum of polarization and its number n are shown in (c) and (d).Due to the SOC, the scaled wrapping interval, [−0.5, 0.5] , is obtained to be half of the interval [−1, 1] .Experimental 85 and theoretical 84 SEPs are included for comparison.
Figure 9. (a) Total, and wrapped partial Berry phases versus .(b) Unwrapped total, and partial Berry phases versus .All the Berry phases are divided by 2π in (a) and (b) so that the wrapping interval is converted from [−π , π] to [− 0.5, 0.5] in (a).(c) Total, wrapped partial, and corresponding spontaneous polarizations versus .The partial polarizations are wrapped into [ −ec/2�, ec/2� ], where ec/� = 59.44 µC/cm 2 is the quantum of polarization.(d) Unwrapped total, partial, and corresponding spontaneous polarizations versus.The unit of polarizations is µC/cm 2 in (c) and (d).Like Fig.6but unlike Fig.4, here, only for simplicity, the auxiliary symbols ∧ and ∼ are not used.All the Berry phases and as a result polarizations are calculated by the standard Berry phase scheme including SP and SOC by GGA+U with U eff = 4 eV for the Multiferroic G-AFM BiFeO 3 along the distortion path as functions of structure from " = 0 " to " = 1 " by step 0.1.The quantum of polarization and its number n are shown in (c) and (d).Due to the SOC, the scaled wrapping interval, [−0.5, 0.5] , is obtained to be half of the interval [−1, 1] .Experimental 85 and theoretical 84 SEPs are included for comparison.

Figure 10 .
Figure 10.Total, wrapped partial, and corresponding spontaneous polarizations versus for (a) LiNbO 3 and (c) LiTaO 3 .The partial polarizations are wrapped into [ −ec/�, ec/� ], where ec/� = 69.80µC/cm 2 is the quantum of polarization.Unwrapped total, partial, and corresponding spontaneous polarizations versus for (b) LiNbO 3 and (d) LiTaO 3 .The unit of polarizations is µC/cm 2 .Like Figs. 6 and 9 but unlike Fig.4, here, only for simplicity, the auxiliary symbols ∧ and ∼ are not used.All the Berry phases and as a result polarizations are calculated by the standard Berry phase scheme including SP by PBE-GGA for the normal FE LiNbO 3 and by LDA for the normal FE LiTaO 3 along the distortion path as functions of structure from " = 0 " to " = 1 " by step 0.1.For comparison, experimental87,87,89 and theoretical84,86,88 SEPs are included in (b) for LiNbO 3 (in (d) for LiTaO 3 ).

Table 2 .
The z components of the Wannier centers, z W n,R , and occupation numbers, n W n,R , as well as summation of n W n,R over the Wannier centers, J (I) =18 n=1 n W n,R , taking the eighteen Wannier centers from n = 1 to J (I) = 18 of class I into account for the NCS R3c (CS R 3 c) hexagonal structure of LiOsO 3 .

of spontaneous polarization of LiOsO 3 by finding the best branch
3.2 of SMs.